% 利用402算例求解网络重构问题的脚本，用来生成样本
% 采用负荷曲线来模拟不同负荷场景的配电网状态
% 考虑分段开关
%% 读取数据
data_read_args.branch_file_path="yantian402/运行方式1-branch_reordered.xlsx";
data_read_args.source_file_path="yantian402/capacity.xlsx";
data_read_args.load_file_path="yantian402/all_load_node_mapped.xlsx";
[buses, branches, substations, init_sw] = yantian_vdata(data_read_args);
Nbr = height(branches);
Nbus = height(buses);

%% 建立图模型
Edgetable = branches(:, {'bus1id', 'bus2id', 'name', 'is_sw'});
Edgetable = mergevars(Edgetable, {'bus1id', 'bus2id'}, "NewVariableName", 'EndNodes');
Edgetable.fixed = ~Edgetable.is_sw;
Edgetable.is_sw = [];
G = graph(Edgetable);

%% 设置初始运行拓扑结构
brid = findedge(G, init_sw.bus1id(~init_sw.status), init_sw.bus2id(~init_sw.status));
initial_grid = true(Nbr, 1);
initial_grid(brid) = false;

%% get load profiles
% 读取负荷曲线并存储为表格
load("yantian402\curve_at_13.mat","needed_curve","sample_ticks");
needed_curve.bus=string(needed_curve.Project_id);
Nt=length(sample_ticks);

%分配负荷曲线
sample_profile_data = struct('P',zeros(Nbus,Nt),'Q',zeros(Nbus,Nt),'model',"real");
for t=1:Nt
    load_at_t=needed_curve(needed_curve.Date==sample_ticks(t),:);
    load_at_t=innerjoin(buses,load_at_t,'Keys','bus','LeftVariables',{'busid'},'RightVariables',{'P','Q'});
    load_at_t=fillmissing(load_at_t,"constant",0,'DataVariables',{'P','Q'});
    sample_profile_data.P(load_at_t.busid,t)=load_at_t.P/1000;
    sample_profile_data.Q(load_at_t.busid,t)=load_at_t.Q/1000;
end

%% 场景1：设置初始拓扑中所有在运支路故障，选择任意支路节点，上优化求解器
% fault_br_name = G.Edges.name(initial_grid);
% fault_br_eid = find(initial_grid);

% 初始化结果存储
% sol = cell(length(fault_br_eid), Nt);
% for k = 1:length(fault_br_eid)
%     this_fault_eid= fault_br_eid(k);
%     parfor m = 1:Nt
%         upbus=buses;
%         % 更新bus负荷并求解
%         if sample_profile_data.model=="real"
%             upbus.Pd = sample_profile_data.P(upbus.busid,m);
%             upbus.Qd = sample_profile_data.Q(upbus.busid,m);
%         end
% 
%         sol{k,m} = reconstruct_branch_fault(G, this_fault_eid, initial_grid, upbus, branches, substations);
%     end
% end

% 存储结果
% save solution_ld_time.mat sol initial_grid G fault_br_name sample_profile_data;

%% 场景2：将相同转供方案的支路放入同一个桶，保证样本总数在桶上均匀分布，也就是说来源于每个桶的样本数量是一样多的
%使用抽样轮数代替时间,总样本数=桶数×抽样轮数
buckets=struct('next',1,'eids',[]);
buckets(1).eids=[305,309];
buckets(2).eids=[101,379];
[buckets.next]=deal(1);
Nbucket=numel(buckets);

Nround=3;
Nround=min(Nround,Nt);
fault_br_eids=zeros(Nbucket,Nround);
sol = cell(Nbucket, Nround);
for k = 1:Nround

    for ii=1:numel(buckets)
        %抽样
        this_fault_eid =buckets(ii).eids(buckets(ii).next);
        fault_br_eids(ii,k) = this_fault_eid;
        buckets(ii).next=buckets(ii).next+1;
        if buckets(ii).next> numel(buckets(ii).eids)
            buckets(ii).next=1;
        end

        %求解
        upbus=buses;

        % 更新bus负荷
        upbus.Pd = sample_profile_data.P(upbus.busid,k);
        upbus.Qd = sample_profile_data.Q(upbus.busid,k);
        
        %发现这个程序暂时没有区分功率倒送情况，再说
        sol{ii,k} = reconstruct_branch_fault(G, this_fault_eid, initial_grid, upbus, branches, substations);
    end

end

% 存储结果
save solution_ld_time.mat sol initial_grid G fault_br_eids sample_profile_data;
